%%%%%%%%%%%%%%%%%% Binomial distribution CI1 %%%%%%%%%%%%%%%%%%%%%%%% 

clear
%format long;
%n=input('Please input N = ')
warning off    % 關掉warning的指令
Result=[];
%for n=5:5:50
n=10;
alpha=0.05;
beta=0.9;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % CI1
 k_item=[0:n];
 p_hat=k_item/n;
 qnorm=norminv(1-alpha/2,0,1);
 p1_tilde=p_hat-qnorm*sqrt(p_hat.*(1-p_hat)/n);
 p2_tilde=p_hat+qnorm*sqrt(p_hat.*(1-p_hat)/n);
 p_tilde=[max(p1_tilde,0);min(p2_tilde,1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%找出I,J=>解出P%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

II=[];
JJ=[];
P=[];
k_item=[0:n];
for k=0:n      
     p1=p_tilde(1,k+1);
     p2=p_tilde(2,k+1);
     
     H1=binocdf(n,n,p1)-binocdf((k_item-1),n,p1)-(1+beta)/2;
     H1(1)=binocdf(n,n,p1)-(1+beta)/2;
     H2=binocdf(k_item,n,p2)-(1+beta)/2; 
     
     c1=min(H1(H1>=0));
     c2=min(H2(H2>=0));
% tolerance interval (c1,c2)
     I=find(H1==c1)-1;
     I=max(I);
     J=find(H2==c2)-1;
     J=min(J);
  
     II=[II,I];    JJ=[JJ,J];
     
%if (I==0 & J==n)
%     p=[NaN NaN];
%else
     bino=[];
     %for x=(I+1):J
     for x=I:J
         strf=[num2str(nchoosek(n,x)) '*p^' num2str(x)  '*(1-p)^' num2str(n-x) '+'];
         bino=[bino,strf];
     end
     bino(end)=[];
     bino=[bino '-' num2str(beta)];
     
%---------------- 
     fun = fcnchk(bino);
     a=feval(fun,0);
     b=feval(fun,1);
     bino(bino=='p')='x';
     
     
     if(a*b<0) %一個解
          p=fzero(bino,[0 1]);
          p=[p NaN];
     else     %無解或兩個解
          for i=0:0.001:1,
              if(feval(fun,i)>0)
                  break;
              end
          end
          if(i==1)|(i==0)
              p=[NaN NaN];
          else
             p1 = fzero(bino,[0 i]);
             p2 = fzero(bino,[i 1]);
             p=[p1 p2];
          end
     end
%end   
     
     
     P=[P;p];
end 
   sprintf('      n  observation  tolerance  tolerance   solution    solution \n                        lower      upper        p           q   \n                        bound      bound  ') 
   output= [n*ones(n+1,1) (0:n)' II' JJ' P]
   sprintf('<Remark>  NaN means no solutions')
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%對解出的P做排序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:length(output(:,1))
   if isnan(output(i,6))==0
      output=[output; output(i,1:4), output(i,6), NaN]; 
   end
end
[B,IX]=sort(output(:,5));
  sprintf('        n         k        I        J       Solve P   ' ) ;
  output=output(IX,1:5);
%output=[n 0 0 0 NaN; output; n n n n NaN];

%%%%%%%%%%%%%%%%%對某個P..找出對應的x, 算出Conf._Level.%%%%%%%%%%%%%%%%%%%%%%%% 
aa=[];
for i=1:length(output(:,1)) %  i=2 ;
   p=output(i,5);
   x=[];
   for j=1:length(output(:,1)) % j=1;
        k=output(j,2);  I=output(j,3); J=output(j,4);
       
      if(I==0)
        temp=sum(binocdf(J,n,p));
      else 
          %temp=sum(binocdf(J,n,p))-sum(binocdf(I,n,p));
        temp=sum(binocdf(J,n,p))-sum(binocdf(I-1,n,p));
      end       
        
      if  temp-beta>10^-3         
          x=[x;k];
      end
   end    
   a=sum(binopdf(min(x):max(x),n,p));
   aa=[aa,a];
   
end
%minconf=min(aa(2:(length(aa)-1))); % 排除k=0 and k=n 沒有解的情況
minconf=min(aa); 
sprintf('      n  observation  tolerance  tolerance   theta    coverage \n                        lower      upper             probability   \n                        bound      bound  ') 
output=[output,aa']

%%%%%%%%%%%%%%%%%%%%%%%%%計算ave. cov.prob.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

P=[];
bd=[];
tempbd=[0;output(:,5);1];
for i=1:length(tempbd) 
   if isnan(tempbd(i))==0
      bd=[bd; tempbd(i)]; 
   end
end

for i=(1:length(bd)-1)
   tempP=(bd(i)+bd(i+1))/2;
   P=[P;tempP];
end

aa=[];
result=[];
for i=1:length(P) %  i=3;
   p=P(i);

   x=[];
   for j=1:length(output(:,1)) % j=1;
        k=output(j,2);  I=output(j,3); J=output(j,4);
        
      if(I==0)
        temp=sum(binocdf(J,n,p));
      else 
          %temp=sum(binocdf(J,n,p))-sum(binocdf(I,n,p));
        temp=sum(binocdf(J,n,p))-sum(binocdf(I-1,n,p));
      end 
      
        if  temp-beta>10^-3;    
            x=[x;k];
        end
   end    
   a=[p,min(x),max(x)];
   
   
   if length(a)==1
       aa=[aa;a NaN NaN];
       resultt=0;
   else 
       aa=[aa;a];
       bino=[];
       for x=a(2):a(3)
              strf=[num2str(nchoosek(n,x)) '*p^' num2str(x)  '*(1-p)^' num2str(n-x) '+'];
              bino=[bino,strf];
       end
       bino(end)=[];  
% average coverage probability for bounded parameter space (bound1, bound2)
      % bound2=0.4;
       %bound1=0.154;
       %newdl=max(bound1,min(bd(i),bound2));
       %newdu=max(bound1, min(bd(i+1),bound2));
       %resultt=1/(bound2-bound1)*int(bino,newdl,newdu);

% average coverage probability for parameter space (0,1)    
  resultt=int( bino,bd(i),bd(i+1));
   end   
   
   result=[result,resultt];    
   
end
%compute the bounds of the function h corresponding to the middle point of two solutions
sprintf('   inter.p     min(x)      max(x)   ' )
aa 
average_coverage_probability =double(sum(result))
result=sprintf('n= %g,  minimum coverage probability =%.4f,  average coverage probability= %.4f',n,minconf,average_coverage_probability)
% result=[n minconf double(sum(result))]
% Result=[Result; result];
% end
% Result